Modelado
Proyecto Final
Alfredo Lefranc Flores, Cynthia Raquel Valdivia Tirado, Rafael Sandoval Fernández, Marco Antonio Ramos Juárez y Francisco Velazquez Guadarrama
2021
Ambiente
En esta segunda parte del proyecto pasaremos a estimar los modelos que mejor puedan predecir la duración de un viaje en bicicleta. Esto con el proposito de poder proveer tanto a la empresa Wheelie Wonka como a usuarios de alguna app de movilidad con información confiable del flujo de bicicletas y bicicletas disponibles para una hora del día determinada.
El primer paso para el análisis es importar los datos procesados de la etapa de EDA. Para ello, cargamos la base y revisamos sus variables. Notamos outliers en la variable de duración. Los truncamos al 99.9% para que no afecten las predicciones de nuestros modelos.
Para poder aumentar la eficiencia de los modelos, transformamos la base a formato ancho. Esto lo hacemos para dos variables con muchas categorías, las estaciones de salida y llegada.
# formato wide de start y end station
stations_wide <- trips_filter %>% select("seq_id", "strt_statn",
"end_statn") %>% as.data.table() %>% one_hot()
trips_filter <- inner_join(trips_filter, stations_wide, by = "seq_id")
trips_filter <- trips_filter %>% select(-c("strt_statn", "end_statn")) %>%
arrange(seq_id)Modelos
Primero separamos la base de validación y de entrenamiento. Usamos una división 80-20. Para esto decidimos estratificar por días y no por viajes. De esta manera tenemos un panorama completo de los días en entrenamiento, y evaluaremos en días nuevos en la base de validación. Esto nos parece un ejercicio más realista dados los objetivos del negocio. De lo contrario, una separación naive correría el riesgo que separarnos los flujos de un día en específico en ambas bases.
Para tener un balance en nuestros sets de entrenamiento y validación, estratificamos por año, mes y día de la semana. Una vez hecha la estratificación, verificamos el balance entre entrenamiento y validación.
# para ello, nos fijamos en todos los días donde tenemos
# viajes
dates <- trips_filter %>% select(c("start_date", "start_month",
"start_weekday", "start_year")) %>% unique()
# estratificamos
treat_id <- treatment_assign(data = dates, share_control = 0.2,
n_t = 1, strata_varlist = c("start_month", "start_weekday",
"start_year"), missfits = "global", seed = 2021, key = "start_date")$data
# pegamos la variable que asigna tratamiento a la base
trips_filter <- left_join(trips_filter, treat_id, by = ("start_date"))
trips_filter <- trips_filter %>% select(-c("strata", "missfit"))Al hacer la separación, nos cercioramos de que exista balance entre ambas bases. Esto lo hacemos con prubas t y F para las variables numéricas.
# summary(trips_filter$treat) el share control se mantiene
# estable
# revisamos balance de clases no estratificadas
# pruebas t en las numéricas
is_num <- lapply(trips_filter, function(x) is.numeric(x))
trips_num <- trips_filter %>% select((names(is_num[is_num ==
TRUE])))
balance_t <- balance_table(data = sample_n(trips_num, 10000),
treatment = "treat")
balance_t$p_value1 <- balance_t$p_value1 %>% round(3)
# filtramos las variables donde se rechaza la hipótesis de
# igualdad al 10%
balance_t_filter <- balance_t %>% filter(p_value1 < 0.1)
datatable(balance_t_filter, options = list(pageLength = 5))De la tabla anterior podemos concluir que 37 variables muestran desbalance, dummies de estaciones y edad. La significancia conjunta no se rechaza. Por ende, consideramos que no es preocupante.
# pruebas F
balance_f <- balance_regression(data = sample_n(trips_num, 100000),
treatment = "treat")
kable(balance_f$F_test, align = "c", booktabs = T, digits = 2,
caption = "Regresión de balance", col.names = c("Estadístico",
"Valor")) %>% kable_styling(position = "center")| Estadístico | Valor |
|---|---|
| F-statistic | 1.19 |
| k | 291.00 |
| n-k-1 | 99708.00 |
| F_critical | 0.87 |
| p_value | 0.01 |
| R cuadrada | 0.00 |
Ahora procedemos a separar las bases en entrenamiento y validación. Acotamos la de entrenamiento para tener una muestra más pequeña. Asimismo, convertimos las bases de covariables en matrices ralas.
train <- trips_filter[trips_filter$treat == 1, ]
y_train <- train$duration[train$treat == 1]
train <- train %>% select(-c("treat", "seq_id", "duration"))
sparse_train <- sparse.model.matrix(~. + 0, data = train)
# validacion
test <- trips_filter[trips_filter$treat == 0, ] %>% select(-c("treat",
"seq_id", "duration"))
sparse_test <- sparse.model.matrix(~. + 0, data = test)
y_test <- trips_filter[trips_filter$treat == 0, ] %>% select(c("seq_id",
"duration"))Una vez realizado el preprocesamiento adecuado procedemos a estimar los modelos: Lasso; Random Forest y XG Boosting.
Lasso
El primer modelo que estimaeros es un LASSO con Cross Validation.
lasso <- cv.gamlr(x = sparse_train, y = y_train, verb = T, nfold = 10)
seg100
5.554548 Guardamos las predicciones del modelo.
Random Forest
El segundo modelo que deseamos probar es el random forest. Probamos con distintas B, desde 100 hasta 700, siendo 500 el óptimo por reducir el error OOS.
# ajuste del set de entrenamiento
dates <- train %>% select(c("start_date", "start_month", "start_weekday",
"start_year")) %>% unique()
# estratificamos
treat_id <- treatment_assign(data = dates, share_control = 0.75,
n_t = 1, strata_varlist = c("start_month", "start_weekday",
"start_year"), missfits = "global", seed = 2021, key = "start_date")$data
# pegamos la variable que asigna tratamiento a la base
train_cf <- left_join(cbind(train, y_train), treat_id, by = ("start_date"))
train_cf <- train_cf %>% select(-c("strata", "missfit"))
# subset de la base para los días elegidos
train_cf <- train_cf[train_cf$treat == 1, ]
y_train_cf <- train_cf$y_train[train_cf$treat == 1]
train_cf <- train_cf %>% select(-c("treat", "y_train"))
sparse_train_cf <- sparse.model.matrix(~. + 0, data = train_cf)Una vez hechos los ajustes, corremos el RF.
# para correr modelos en paralelo
cores <- detectCores()
cl <- makeCluster(cores)
inicio <- Sys.time()
# c(100,200,350,500,700)
random_forest <- map(c(500), function(z) ranger(y = y_train_cf,
x = sparse_train_cf, num.trees = z, mtry = ncol(train) %>%
sqrt() %>% floor(), min.node.size = 1, importance = "impurity",
status.variable.name = 1))
(tiempo <- Sys.time() - inicio)
stopCluster(cl)
error_prediccion <- tibble(trees = c(500), oob_error = map_dbl(random_forest,
~.$prediction.error))
# error medido en horas ggplot(error_prediccion, aes(trees,
# oob_error/60)) + geom_point() + geom_path() + theme_bw()df_imp <- data.frame(names = random_forest[[1]][["variable.importance"]] %>%
names(), importance = random_forest[[1]][["variable.importance"]]) %>%
arrange(-importance)
# predicciones OOS
y_test$pred_rf <- predict(random_forest[[1]], data = sparse_test)$predictionsPredicting.. Progress: 48%. Estimated remaining time: 33 seconds.y_test <- y_test %>% mutate(oos_rf = (duration - pred_rf)^2)XGB
Finalmente, probamos un XGB. Entrenamos el modelo y unimos las predicciones a la base de evaluación.
# Preparar la base de entrenamiento y de validación
dtrain <- xgb.DMatrix(sparse_train, label = y_train)
dtest <- xgb.DMatrix(sparse_test, label = y_test$duration)
watchlist <- list(train = dtrain, eval = dtest)
# Entrenamiento del modelo
param <- list(max_depth = 5, learning_rate = 0.06, objective = "reg:squarederror",
eval_metric = "rmse", subsample = 0.85, colsample_bytree = 0.7)
xgb_model <- xgb.train(params = param, dtrain, early_stopping_rounds = 10,
nrounds = 100, watchlist)# Predicción
y_test$pred_xgb <- predict(xgb_model, sparse_test)
y_test <- y_test %>% mutate(oos_xgb = (duration - pred_xgb)^2)Una vez estimados los tres modelos, procedemos a evaluar su desempeño, tanto para métricas tradicionales como para métricas más orientadas a los fines del negocio.
Desempeño de modelos
Medidas tradicionales
Como medidas tradicionales, al ser la variable objetico numérica (la predicción de la duración de un viaje en particular), podemos emplear el R^2 y además podemos comparar ambos modelos midiendo el tamaño de los residuales
| \[R^2\] | Minutos residuales totales | Desviación en minutos | |
|---|---|---|---|
| RF | 0.116 | 22013847327 | 0.179 |
| Lasso | 0.073 | 23082505149 | 0.344 |
| XGB | 0.134 | 21568657965 | 0.142 |
Medidas alternativas: aplicación de negocio
Predicciones y flujo de biciletas por estación
Para mostrar el desempeño de los modelos, modelaremos los flujos de bicicleta de los días que componen nuestra base de validación pero ahora en lugar de usar la duración del viaje real, usaremos la predicción de cada modelo. Para ellos crearemos primero una base de datos que compare los flujos de viaje cada 10 minutos entre lo que en realidad pasó y lo que predice nuestro modelo.
En la siguiente tabla se muestran los flujos acumulados de bicicletas que salen y llegan (con datos reales y predecidos en el set de validación). Por motivos de costo computacional, solamente mostramos algunos días completos para algunas estaciones aleatorias. Por favor haga click primer primero en estación y luego en fecha para que la tabla esté ordenada de acuerdo a la hora y de esta manera se pueda apreciar la acumulación de viajes.
Predicción de flujos de llegadas de bicicleta para cada modelo
# Deployment de las bases de datos en una app (selecciono de
# manera aleatoria algunas debido al peso computacional)
random_stations <- sample(3:133, 5, replace = FALSE)
output <- plantilla_llena %>% filter(strt_statn %in% random_stations)
# output<-sample_n(plantilla_llena, 1000)
output <- output %>% select(date, hour, strt_statn, acum_salida,
acum_llegada_real, acum_llegada_pred_rf, acum_llegada_pred_lasso,
acum_llegada_pred_xgb)
names(output) <- c("Fecha", "Hora", "Estación", "Salidas", "Llegadas(real)",
"RF", "Lasso", "XGB")
datatable(output, options = list(pageLength = 15))Error promedio a la hora de corte
Derivado de estas tablas, se nos ocurrió que una forma muy util de probar el desempeño de los modelos es a través del error del flujo en determinado punto del día. Esto pensando en dos escencarios: primero, que en las operaciones de este tipo, hay un punto de corte en el que las estaciones se tienen que llenar o vaciar de bicicletas; segundo, que tal vez como usuario, el saber cuantas bicis hay disponibles en cada estación es muy útil. Para este ejemplo, propusimos un punto de corte (de rellenado o vaciado de bicicletas) a las 12 de la noche, pero bien podría ser en la madrugada, o en la tarde. En la siguiente tabla se muestra el error promedio en el cálculo del flujo de bicis para la hora de corte.
Predicción de bicicletas por estación con otro acercamiento
Una vez cubiertos los primeros dos objetivos del trabajo, tuvimos que decidir cómo resolver el tercero: hacer una predicción en ventanas de 10 minutos de cuántas bicis habría en cada estación. Para hacer esto con nuestra predicción de duraciones, se nos ocurrían dos opciones.
Simular datos de días no observados para, en tiempo real y dado un número de bicicletas en cada estación y los viajes activos, predecir cuáles terminarían su viaje en ese lapso así como el destino. Si bien esta opción nos parecía útil desde el punto de vista de las necesidades del negocio, creímos que implicaba otra naturaleza de análisis independiente del realizado hasta ahora.
DESCRIBIR LA OPCIÓN TOMADA
En esta sección se presenta un modelo distinto para predecir la disponibilidad de bicicletas por fecha y hora. Para ello, se aplica una transformación a los datos de manera que estos reflejen el flujo de entradas y salidas por cada intervalo de diez minutos. La idea es que la predicción de dichos flujos a futuro pueden ayudar a notar, por ejemplo, en qué estación han salido más bicicletas de las que han entrado, indicando una baja disponibilidad. Para ello, las observaciones se agregan por coincidencia de fecha, día, hora y estación. Naturalmente, una estación que no registra un flujo no tuvo entradas ni salidas de bicicletas para ese grupo horario. Con esto obtenemos una base de datos reales que cuenta con una observación para cada intervalo del día, y para todos los días de los que se tiene registro. Con esto se intenta predecir los flujos a futuro posteriormente.
# Reestructuración de la base de datos
load("flux_df.RData")
flux_df <- flux_df %>% select(c("start_date", "end_date", "strt_statn",
"start_month", "start_weekday", "start_hour", "start_minute",
"end_statn", "end_month", "end_weekday", "end_hour", "end_minute"))
# Se transforma a los minutos por intervalos de 10 minutos
flux_df$start_minute <- cut(flux_df$start_minute, breaks = seq(0,
60, 10), labels = seq(1:6), right = F)
flux_df$end_minute <- cut(flux_df$end_minute, breaks = seq(0,
60, 10), labels = seq(1:6), right = F)
# Columnas auxiliares para medir el flujo
flux_df$arrival <- 1
flux_df$departure <- 1
# Se agrega todo movimiento de entrada y salida por estación,
# fecha y hora
outflux <- aggregate(departure ~ start_date + strt_statn + start_hour +
start_minute, FUN = sum, drop = F, data = flux_df)
influx <- aggregate(arrival ~ end_date + end_statn + end_hour +
end_minute, FUN = sum, drop = F, data = flux_df)
rm(flux_df)
outflux[is.na(outflux$departure), "departure"] <- 0
influx[is.na(influx$arrival), "arrival"] <- 0
colnames(outflux) <- c("date", "station", "hour", "minute", "outflux")
colnames(influx) <- c("date", "station", "hour", "minute", "influx")
# Indicadoras para la unión
outflux$ID <- with(outflux, paste0(date, " ", station, "-", hour,
"-", minute))
influx$ID <- with(influx, paste0(date, " ", station, "-", hour,
"-", minute))
# Base con los movimientos de entrada y salida para cada
# grupo horario de cada día
flux <- merge(outflux, influx[, c("ID", "influx")], by = "ID",
all.x = T, all.y = F)
rm(outflux, influx)
# Resultado neto
flux$netRes <- flux$influx - flux$outflux
flux$weekday <- weekdays(flux$date)
flux <- flux %>% select(-c("date", "outflux", "influx"))
# Preparación de variables
flux$station <- flux$station %>% as.factor
flux$hour <- flux$hour %>% as.factor
flux$minute <- flux$minute %>% as.factor
flux$netRes <- flux$netRes %>% as.numeric
flux$weekday <- flux$weekday %>% as.factor
# Se guarda la base save(flux, file='flux.RData')# Separación en test-train, y entrenamiento del modelo
#load('flux.RData') # Para cargar poder correr la generación de datos por separado.
flux$ID <- NULL
set.seed(2020)
sample <- sample.split(flux$weekday, SplitRatio = .75)
train_f <- subset(flux, sample == TRUE) # Train
test_f <- subset(flux, sample == FALSE) # Test
test_target <- select(test_f, 'netRes')
test_f$netRes <- NULL
# Modelo: Random Forest
cores <- detectCores() # Paralelización
cl <- makeCluster(cores)
inicio <- Sys.time()
flux_forest <- ranger(netRes ~ .,
#data = train_f,
data = under,
num.trees = 300,
classification = T,
verbose = T,
oob.error = F,
write.forest = T)
tiempo <- Sys.time() - inicio
stopCluster(cl)
# Se guarda el modelo
save(flux_forest, file = "flux_forest.RData")
# Evaluación
flux_preds <- predict(flux_forest,
data = test_f,
type = "response",
verbose = T)
comparativa <- data.frame(test_target, flux_preds$predictions)Finalmente, los resultados de este último modelo se muestras a continuación:
Conclusiones
Con estos resultados es posible optimizar la gestión del servicio de bicicletas otorgado por la empresa Wheelie Wonka. Este pronóstico preciso sirve de guía para que los pasajeros puedan organizar mejor su hora de salida y las rutas de desplazamiento. Además, es beneficioso para el proveedor de servicios de bicicletas compartidas Wheelie Wonka en términos de mejorar la satisfacción de sus clientes y organizar con efectividad el horario de entrega y distribución de bicicletas y así eficientar sus recursos y disminuir pérdidas.
De manera más técnica, el modelo que resultó dar los mejores resultados tanto desde una perspectiva de análisis de residuales tradicional como desde un problema aplicado de negocios fue el XGB.